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Abstract 

We investigate the expansion dynamics of a dilute Fermi gas at unitarity in the context of 
dissipative fluid dynamics. Our aim is to quantify the effects of shear viscosity on the time evolution 
of the system. We compare exact numerical solutions of the equations of viscous hydrodynamics to 
various approximations that have been proposed in the literature. Our main findings are: i) Shear 
viscosity leads to characteristic features in the expansion dynamics; ii) a quantitative description 
of these effects has to include reheating; iii) dissipative effects are not sensitive to the equation 
of state P(n, T) as long as the universal relation P = |£ is satisfied; iv) the expansion dynamics 
mainly constrains the cloud average of the shear viscosity. 
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I. INTRODUCTION 



A cold, dilute Fermi gas in which the scattering length can be tuned to infinity by means 
of a Feshbach resonance provides a new realization of a strongly correlated quantum fluid []]- 
[J. On resonance the two-body scattering amplitude saturates the s-wave unitarity bound, 
and the corresponding many-body system is referred to as the Fermi gas at unitarity. An im- 
portant manifestation of strong correlations is the observation of nearly ideal hydrodynamic 
flow [4|. In the Fermi gas at unitarity nearly ideal flow was observed in samples containing 
as few as 10 5 atoms, with an interparticle spacing between the atoms on the order of several 
10 3 A, much larger than the range of the interaction. 

This result implies that dissipative effects must be very small In a normal fluid 
dissipative phenomena are governed by three transport coefficients, the shear viscosity rj, 
the bulk viscosity (, and the thermal conductivity k. The Fermi gas at unitarity is scale 
invariant and the bulk viscosity is zero 

aa. 

Most of the experiments conducted so far, 
like collective oscillations and the expansion from a deformed trap, involve scaling flows. 
In these flows the cloud remains nearly isothermal and the experiment is not sensitive to 
thermal conductivity js, 9]. The observation of nearly ideal flow therefore requires that the 
shear viscosity is very small. 

From a theoretical point of view we know that the shear viscosity of a weakly interacting 
gas scales as rj ~ p/c, where p is the mean momentum and a is the scattering cross section. 
This implies that the shear viscosity of a strongly interacting gas is expected to be small. 
At unitarity the cross section is a = 4ir/k 2 , where k is the momentum transfer. The average 
cross section in a thermal gas is a ~ 4^mT, and the shear viscosity is 77 ~ (mT) 3 / 2 
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111 ] . This result is reliable as long as T is much larger than the critical temperature T c 



for superfluidity. Below T c the nature of the excitations changes and rj ~ 1/T 5 [12]. These 
results indicate that the shear viscosity has a minimum at a temperature on the order of 
T c , but kinetic theory cannot reliably predict how small the minimum value of the shear 
viscosity is. The region near T c can be studied using sum rules jl3j ]. or with the help Kubo's 
formula and many body perturbation theory [jjj]. 

It has been argued that the quantum mechanical uncertainty relation implies a lower 
bound 77/77, ~ pl m fp > h [151]. Here, n is the density, l m f p ~ l/{na) is the mean free path, and 
h is Planck's constant. A more precise bound has emerged from the study of holographic 
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dualities in string theory. Kovtun, Son, and Starinets proposed that the shear viscosity to 
entropy density ratio rj/s is bounded from below by ^/(47t/cb) 16|, where ks is Boltzmann's 
constant. 

Simple estimates of the shear viscosity to entropy density ratio based on experimental 
data indicate that 77/s in the unitary Fermi gas is indeed close to the proposed bound. The 
:irst attempt to determine rj/s from data was based on the damping of collective modes 



17H19|. The damping constant can be related to the ratio E/E, where E is the total energy 



of the mode and E is the energy dissipated by viscous effects. For a scaling flow E is 
proportional to the spatial integral of the shear viscosity. It was found that the ratio of the 



trap averages of 77 and s has a minimum value of (rj)/ (s) ~ 0.5 [18|, [l9|, where we have set 
U — ks — 1. 

More recent analyses are based on the dynamics of an expanding cloud Q, Q-22]. These 
studies utilize approximate solutions of the equations of dissipative hydrodynamics. The 
first approximations is that entropy is assumed to be conserved, which is equivalent to the 
assumption that the system is in contact with an energy sink that removes the heat gen- 
erated by dissipative effects. The second approximation is that the Navier-Stokes equation 
is converted into a set of ordinary differential equations by taking moments. If only the 
lowest moments are included then the result is only sensitive to the spatial integral of the 
shear viscosity. An analysis of the expansion of a rotating cloud gives values as small as 
(v) I i s ) — 0-2 22[- A more refined approximation was used in |23| to analyze the high 
temperature limit of the shear viscosity. We will describe this method in Sect. |VJ 

The goal of our present work is to test these approximations by performing numerical 
studies of the equations of viscous hydrodynamics for an expanding cloud of a dilute Fermi 
gas at unitarity. This paper is structured as follows. In Sect. HT1 and ITTTl we introduce the 
equations of dissipative hydrodynamics for a scale invariant non-relativistic fluid. We discuss 
exact solutions for the ideal case and approximate solutions for the viscous case in Sects. [TV] 
and|Vl Numerical solutions are discussed in Sect. |VT]and our conclusions are summarized 
in Sect. EID 
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II. DISSIPATIVE HYDRODYNAMICS 



We will consider the unitary Fermi gas in the normal phase. In that case there are five 
hydrodynamic variables, the mass density p, the flow velocity v, and the energy density S. 
These variables satisfy four hydrodynamic equations, the continuity equation, the Navier- 
Stokes equation, and the equation of energy conservation, 

% + V • (pv) = 0, (1) 

^ + v,n tf = o, (2) 

§ + v.r = o. (3) 

The total energy density is the sum of the internal energy density and kinetic energy density, 
S = £ + \pv 2 . These equations close once we supply constitutive relations for the stress 
tensor H^- and the energy current j/ as well as an equation of state. The unitary Fermi gas 
is scale invariant and the equation of state is P = |£o- The stress tensor is given by 

Uij = pViVj + P5ij + 5U i:j , (4) 

where SUij is the dissipative part. The dissipative contribution to the stress tensor is SUy = 
—r\Gij with 

Oij = (ViVj + VjVi - -5ij{V k v k )^ , (5) 

where r\ is the shear viscosity and we have used the fact that the bulk viscosity of the unitary 
Fermi gas is zero. The energy current is 

Ji = ViW + 5j t e , (6) 

where w = So + P is the enthalpy density. The dissipative energy current is 

Sj^SIUjVj-KViT, (7) 

where T is the temperature and k is the thermal conductivity. We note that the temper- 
ature T = T(n,P) is a function of the density n = p/m and the pressure. In order to 
determine T we need the equation of state in the form P = P(n,T). Universality implies 
that P(n,T) = m~ 1 n 5 / 3 f n (mT/n' 2 / 3 ) where f n (x) is a universal function that has to be de- 
termined experimentally or using quantum Monte Carlo methods. The situation simplifies 



in the high temperature limit where P = nT. We provide a parameterization of P(n,T) 
at all temperatures in Appendix |X] Universality also restrict the dependence of the shear 
viscosity and thermal conductivity on the density and the temperature. We can write 

/ TflT \ { TlTlT \ Tl 

v (n,T)=a n ^j n, K {n, T) = a n J - , (8) 

where a n (y) and a n (y) are universal functions of y = mT/n 2 ^ 3 . The relative importance 
of thermal and momentum diffusion can be characterized in terms of a dimensionless ratio 
known an the Prandtl number, Pr = c p t]/(pk), where c p is the specific heat at constant 
pressure. In the high temperature limit c p = p/m and Pr = a n /a n . Kinetic theory predicts 
that in this limit Pr = 2/3 fl. 



III. INITIAL CONDITIONS AND CHOICE OF UNITS 

We will consider the expansion of a dilute Fermi gas after release from a harmonic trap. 
The density distribution in the initial state satisfies the equation of hydrostatic equilibrium, 
VP = — nVU, where V(x) = ^mujfx 2 is the trapping potential. If the gas is isothermal 
then this equation is solved by the local density approximation 

n (x) = n(p(x), T) , n(x) = p — V(x) , (9) 

where n(p, T) is the density in thermal equilibrium. The function n(p, T) can be determined 
from the equation of state as explained in appendix |A] The simplest case is the high 
temperature limit. In this limit the initial density is a Gaussian 

n (x) = n exp > ( 10 ) 

with R 2 = (2T) / (mui 2 ) and n = N/ (n 3 ^ 2 R x R y R z ). The total number of particles is denoted 
by N. In the following we will use a dimensionless coordinate variable X{ = Xi/xo where 

x 2 = — ( 3N 

3m \u x tu y u z 

This variable can be used for any initial condition, not just the Gaussian initial condition 
given in equ. fflO|) . We will focus on axially symmetric traps with u x = u y = u± and 
u z = \u±. In the high temperature limit the dimensionless density n = nxjj is given by 




where Eq is the total (potential and internal) energy of the trapped gas and Ep = 
(WXy^Nux. The central density is n = XN(E F / '£ ) 3/2 /tt 3 / 2 . 

We can write the equations of hydrodynamics in dimensionless variables by introducing 
a scaled time variable t = u±t as well as a scaled velocity, energy density, and pressure, 

£ = ^S, P -%P. (13) 



The scaled mass density is p = px\jm. Using these variables the equations of fluid dynamics, 
equ. (jTJ|3]) , remain unchanged except for the change from dimensionful to dimensionless 
hydrodynamic variables. The dimensionless shear viscosity is fj = XqT]/ \muj±). We can write 
fj = a n n with 

- - 3 1 

Gn ~ 2(3AiV)V3 Q; "' ^ 
where a n is the universal function introduced in equ. (jHJ). Finally, we can introduce a 
dimensionless temperature and chemical potential, T = T/(miu]_xl) and ft = ///(mw|ig). 

IV. EXACT SOLUTIONS 

In order to test the accuracy of the numerical hydrodynamics code we have studied 
a number of exactly solvable test cases. In ideal hydrodynamics there are exact scaling 
solutions for the expansion from rotating and non-rotating traps. Consider a density profile 
of the form 

1 fx 2 y 2 A¥ \ 

n[X,t) ~ b x (t)b y (t)b z (t) [bl(t) + &g(t) + bl(t) ) ' ( 5) 

where F(x) is an arbitrary function and the scale parameters bi(t) satisfy the initial condi- 
tion foj(O) = 1. This ansatz satisfies the continuity equation with a velocity field given by 
Vi(x,t) = ai(t)xi with aij = fej/ftj. The initial condition for the pressure can be determined 
by integrating the equation of hydrostatic equilibrium 

P {x) = - J n {x)VV{x) ■ dx. (16) 

In the limit T 3> Tp the function F(x) is a Gaussian and the initial pressure is determined 
by the ideal gas equation of state, P = nT. In the initial state the temperature is constant 
and the chemical potential is parabolic. In ideal hydrodynamics the evolution of the system 
preserves these properties. The Gibbs-Duhem relation dP = ndp + sdT implies that the 
force (VP)/n = V/x is exactly linear at all times. 



The Euler equation is equivalent to three coupled ordinary differential equations for the 
scale parameters b{. We get jsl, Q] 



oo? 1 



(17) 



(b x b y b z )W k ' 

with the initial conditions 6^(0) = 1 and foj(0) = 0. In the case of axial symmetry this set of 
equations reduces to two independent equations for b± = b x = b y and b z . These differential 
equations have to be solved numerically. 

This solution can be generalized to an initial velocity field that corresponds to a rotating 
trap. The velocity field can be chosen to be irrotational, v = aV(xz), or rigidly rotating, 
v = Qy x x. Here we have chosen the direction of the angular momentum to be in the 
y-direction. In the rotating case the profile function in equ. (IT51) has to be generalized 
to include an off-diagonal xz-term. In total one has to solve for ten functions, the four 
scale parameters b x ,b y , b z and b xz , the chemical potential at the center of the trap, and five 
functions characterizing the velocity field, a x , a y , a z , a and fl The equations of motion are 



given m 
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There are no known exact solutions for an expanding gas in the dissipative case. However, 
we can find scaling solutions to the continuity and Navier-Stokes equation if the local shear 
viscosity is of the form rj = i]oP/T, where 770 is a constant. Note that at high temperature 
P = nT and 77 = rjon. Also, for any scaling solution T = const and Vr? = r]o(VP)/T, which 
implies that both the ideal and the dissipative forces are proportional to the gradient of the 
pressure. The Navier-Stokes equation then leads to the coupled set of equations 

h a;j 2pu± (b± b z \ 

x (bib z y/% ± b ± \b ± bj 

* (bib z y/% z + b z \b ± b z )> [U) 

where we have specialized the solution to the case of axial symmetry and we have defined 
(3 = t]qU±/(3T ), where T is the initial temperature. We can write the parameter (3 as 

a _ (gn) 1 _ 2 (a w ) , , 

' (3N\y/3(E /E F ) 3(E /E F ) } 1 ] 



where (a n ) is the trap average of a n 

1 

-j „ o ^ 2/3 



1 f o / mT \ 
a n) = ir T d x a n \ — -r-wh n oi x ) > (21) 



and we have used equ. (111]) . The solution of equ. (I18II19P does not conserve energy. Instead, 
it satisfies a modified energy equation 

f + V-r = -|K) 2 , (22) 

which contains a sink that removes the heat generated by dissipative effects. This means that 
the scaling solution conserves entropy even if the shear viscosity is not zero. The produced 
entropy is removed by the heat sink. We also note that equ. ( jT8¥T9T) can be generalized to 
the rotating case, see Q]. 



V. VISCOUS HYDRODYNAMICS: APPROXIMATE SOLUTIONS 

If dissipative effects are mostly governed by viscous forces, and reheating is not important, 
then the scaling solution introduced in the previous section is a useful approximation to the 
full hydrodynamic equations. We will show below that for an expanding system this is not 
the case. A more useful approximation was recently proposed in [23J]. We will assume that 
the local shear viscosity is proportional to the density, i] = a n n, where a n is a constant. 
The basic idea is to focus on the force f\ = (VjP)/n rather than the pressure itself. The 
Navier-Stokes equation is 

m g + ,.v) B( = / ( + ^pl. (23) 
With the help of the Navier-Stokes equation the energy equation can be written as 

(| + v • v + \ (v • efj h + (v^u- - ~ (v»v^) £ = -jj^i (24) 

where q = ^((Jij) 2 is the heating rate. The basic idea is to assume that even in the dissipative 
case the velocity field and the force remain linear in the coordinates. If the velocity is linear 
and 7] ~ n then all terms in equ. ( 1231) are linear in Xj. Also, equ. ( 124"]) is independent of 
the pressure and all the remaining terms are linear in Xj. We write /, = ciiXi, Vi = otiXi (no 
sum over i) and use the scaling ansatz ( IT5]) for the density. The continuity equation requires 
a>i = bi/bi. The scale parameters and 6j are determined by the coupled equations 

(25) 
(26) 
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b z 
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b z 



a - = -3°H 5 c + sJ + ^bT\vr^) • (27) 
a » = -3 a 4 4 ^ +2 ^j + ^rW^j • (28) 

where (3 is defined in equ. (120]) . The initial conditions are &j_(0) = & 2 (0) = 1, fej_(0) = & z (0) = 
as before, and a±(0) = uj_, a 2 (0) = We note that the solutions of equ. f l2"5TT2"gj) provide 
an exact solution of the continuity, Navier-Stokes, and energy conservation equation. The 
solution is approximate in the sense that one cannot in general find a consistent expression 
for the pressure P(/i,T) such that P,S,n are related by thermodynamic identities. 



VI. VISCOUS HYDRODYNAMICS: NUMERICAL RESULTS 
A. Numerical tests 

Ideal hydrodynamic simulations were carried out using the VH1 code written by Blondin 
and Lufkin [26]. VH1 uses the PPMLR (Piecewise-Parabolic Method, Lagrangian- Remap) 



method developed by Colella and Woodward 271, |28[. The hydrodynamic equations are 
written in the form of conservation laws and solved in Lagrangian coordinates. A Lagrangian 
time step is followed by a piecewise parabolic remap onto an Eulerian grid. We have modified 
VH1 to include viscous corrections to the stress tensor and the energy current. In the current 
work we have not included the effect of thermal conductivity. In the ideal case the cloud 
remains isothermal during the expansion and VT = 0. Dissipative effects lead to non-zero 
temperature gradients, but if the viscosity and thermal conductivity are small then the 
corresponding correction to the energy current is second order in small quantities. We will 
verify this statement in Sect. IVI Bl 

We first consider the evolution of a Gaussian density profile in ideal hydrodynamics. The 
initial condition is given by equ. (Tl2|) where we have chosen E /Ep = 1 and n = 1 (the 
ideal evolution is independent of the number of particles). The aspect ratio of the cloud is 
A = 0.045. In Fig. [TJwe show the evolution of the density and the velocity. The points are 
the results of the numerical solution using VH1 and the lines are semi-analytic results based 
on solving the coupled set of ordinary differential equations fTTTT) . The numerical calculation 
was performed on a fairly coarse grid with a grid spacing Ax = 0.2. We observe that the 
numerical calculation is nevertheless very accurate. 




FIG. 1: Density and velocity profile of an expanding gas cloud in ideal hydrodynamics. The 
initial condition is a Gaussian density profile with Eq/Ef = 1. The top panel shows the density 
n(x, 0, 0, t) for i = 0.0, 0.25, 0.50, . . . , 1.25. The solid lines show the analytic solution of the Euler 
equation, and the dots show numerical results. The bottom panel shows the velocity v x (x,0,0,t) 
for i = 0.0,0.25,0.50,0.75. 

The lower panel of Fig. [T] shows that the velocity field tracks the linear behavior of the 
analytic solution only up to some maximum distance which slowly grows with time. The 
turnover of the velocity field is related to the fact that we impose a minimum density and 
pressure (10~ 15 of the initial central density and pressure). Beyond the point at which the 
minimum pressure is reached there are no pressure gradients and therefore no acceleration. 
The sharp discontinuity in the velocity does not lead to numerical problems. A smooth 
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FIG. 2: Energy of the expanding gas cloud as a function of time. We show the kinetic, internal, 
and total energy. The solid lines show the analytical result for the ideal evolution, the dashed lines 
show the analytical result without reheating for f3 = 0.066, and the points are from a numerical 
calculation with Eq/Ef = 1 and f\ = a n n with a n = 0.1. In the simulation we remove the heat 
generated by dissipative effects. 

turnover of the velocity field can be achieved by considering slightly modified initial condi- 
tions. If we solve the equation of hydrostatic equilibrium in a potential which is harmonic 
at short distances, but grows as V ~ \x\ a with a < 1 at large distances, then the velocity 
field of the expanding cloud will go to zero smoothly as \x\ — > oo. 

We next consider the dissipative evolution of a Gaussian density profile. We take the 
shear viscosity to be of the form 1] = a n n with a constant a n . In order to compare with the 
solution discussed in Sect. HV]we include a sink in the equation of energy conservation as 
defined in equ. (1221) . Note that we can write the divergence of the dissipative energy current 
as 

V 4 (5j, e ) = V 4 (vjSUij) = ~\ K) 2 + ViVjdUij . (29) 

The first terms corresponds to viscous heating and the second term describes the work done 
by viscous forces. Adding a heat sink implies that we only keep the effect of the work term. 

The evolution of the kinetic, potential, and total energy for a n = 0.1 and Eq/Ef = 1 
are shown in Fig. [2j The solid lines show the ideal evolution determined by equ. (|T71) 
and the dashed lines show the dissipative result given by the solution of equ. (I18|I19I) for 
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FIG. 3: Energy of the expanding gas cloud as a function of time. We show the kinetic, internal, 
and total energy. The solid lines show the analytical result for the ideal evolution, and the dashed 
lines show the analytical result without reheating for /3 = 0.066. The data points show the result 
of a numerical calculation with Eq/Ef = 1 and fj = a n n with a n = 0.1. The simulation includes 
all dissipative terms in the energy current. 

= \ocn{Ep / Eq) = 0.066. The data points come from a numerical calculation based on the 
dissipative version of VH1. We observe that the data agree very well with the analytical 
solution. The main difference between the dissipative and the ideal solution is that a fraction 
of the kinetic energy is converted to heat. The heat is absorbed by the sink and lost to the 
system. The evolution of the internal energy is only affected indirectly, via the effect of 
dissipation on the evolution of the radius of the system. 

B. Effects of reheating 

The evolution of the energy in a complete simulation, including the effects of reheating, 
is shown in Fig. El The system parameters are the same as in the previous section. We 
observe that the total energy is conserved to a very good accuracy. Reheating increases the 
internal energy as compared to the result in ideal hydrodynamics. Hydrodynamic evolution 
converts the added internal energy into kinetic energy. This implies that reheating leads to 
reacceleration. 
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FIG. 4: The top panel shows the transverse and longitudinal scale factors of the expanding gas 
cloud as a function of time. The initial density has a Gaussian profile with A = 0.05 and Eq/Ef = 1. 
The viscosity is of the form fj = a n n with a n = 0.1. The solid lines show the analytical result for 
the ideal evolution, the dashed lines correspond to the dissipative solution without reheating, and 
the dotted line is the approximate solution described in Sect. |V] The points are from a numerical 
calculation. The bottom panel shows the ratio of the transverse scale factor over the result in ideal 
hydrodynamics for the calculations shown in the top panel. 

This is shown in more detail in Fig. HI The upper panel shows the time evolution of 
the Gaussian radii R±(t) and R z (t). We have normalized R±(0) = 1 and R z (0) = 1/A so 
that R±(T) = R z (t) corresponds to an aspect ratio of one. The solid lines show the result 
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FIG. 5: This figure shows the time evolution of the aspect ratio Ar = \R±(t)/R z (t). The solid 
lines show the analytical result for the ideal evolution, the dashed lines correspond to the dissipative 
solution without reheating, and the dotted line is the approximate solution described in Sect. [VJ 
The points are from a numerical calculation. All curves correspond to a Gaussian initial condition, 
and the viscosity is of the form fj = a n n with a n = 0.1. 

in ideal hydrodynamics, the dashed lines show the result without reheating, the dotted line 
shows the approximate solution including reheating discussed in Sect. [V] and the data points 
are obtained from a numerical calculation. The lower panel shows the ratio R±{t)/R l f{t) 
where R l f(t) is the transverse size in ideal hydrodynamics. We observe that the two radii 
initially track the prediction of the calculation without reheating: Viscosity slows down the 
expansion in the short direction, and accelerates the system in the longitudinal direction. At 
later times reheating leads to an acceleration in both directions. The lower panel of Fig. H] 
shows that the transverse size almost goes back to the prediction of ideal hydrodynamics. 
For a shear viscosity which is linear in the density this behavior is very well described by 
the linear force model discussed in Sect. [V] 

Fig. |5] shows that even if reheating is taken into account shear viscosity leads to significant 
effects in the evolution of the aspect ratio Ar = XR±/R Z . The shape of A R (t) is similar 
in the model without reheating and the numerical simulation, but the magnitude of the 
dissipative effect is about a factor 2 smaller if reheating is taken into account. We observe 
that shear viscosity leads to a characteristic bending of A R (t) in the regime A R ~ 1. In ideal 
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FIG. 6: This figure shows the viscous correction to the crossing time t cr defined by An(t cr ) = 1. 
The top panel shows 5i cr as a function of the viscosity a n , where f) = a n n, for E^/Ep = 1 and a 
fixed initial aspect ratio A = 0.05. The solid line is the result of the dissipative calculation without 
reheating and the points are from a numerical calculation. The bottom panel shows St cr as a 
function of A for a fixed viscosity a n = 0.1 

hydrodynamics acceleration takes place at early times i < 3. Dissipative forces and reheating 
lead to longitudinal acceleration which occurs on a much longer time scale. Observing this 
behavior not only constrains the value of the shear viscosity, it also demonstrates that the 
systems continues to behave hydrodynamically even at very late times. 

The effect of shear viscosity on the evolution of Ar(£) can be quantified in terms of the 
"crossing time" t cr defined by An(t cr ) = 1. Viscosity leads to a shift St cr in the crossing 
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FIG. 7: This figure shows the evolution of the temperature profile in viscous hydrodynamics. The 
data points show the temperature T(x,0,0,t) determined in a numerical simulation with a = 0.1 
at several different times i = 0, 1.68, 2.70. The lines are the result in ideal hydrodynamics. 

time as compared to ideal hydrodynamics. In Fig. we show 5i cr as a function of a n and 



in the limit A 1 up to higher order corrections in A. Neglecting the effects of reheating 
the correction to the crossing time is jg] 



where (a n ) is the average of a n over the initial density of the trap. The change in t cr if 
reheating is included is shown in Fig. EJ In the upper panel we show the dependence of 
5i cr on a n . We observe that the effect remains linear if reheating is included, but that the 
sensitivity of 5i cr to a n is reduced by about a factor of 2. The lower panel shows that 
the correction factor depends on the geometry. If reheating is included then the sensitivity 
of 5i cr to the shear viscosity becomes very small in the limit of strongly deformed traps 
(A — > 0). This is related to the fact that in this limit all internal energy is converted to 
transverse motion, irrespective of whether or not there is dissipation. 

Fig. [7] shows the effect of reheating on the temperature profile of the cloud. The solid 
line shows the temperature profile at different times during the ideal evolution, and the 
data points come from a numerical simulation with a n = 0.1. The increase in the average 



A. Ideal hydrodynamics predicts that t 



y/l/(^u±) with 7 = 2/3. This result is correct 




(30) 
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FIG. 8: This figure shows the evolution of the produced entropy per particle AS/N as a function 
of time. The thick line shows the result in a numerical simulation with Eq/Ef = 1 and a n = 0.1 
and the thine line is the entropy removed by the heat sink in the model discussed in Sect. IIV1 

temperature relative to the result in ideal hydrodynamics is due to reheating combined with 
a decrease in the expansion rate. In our simulation we have used a spatially constant a n 
and the high temperature equation of state P = nT. In this case both the dissipated heat 
and the specific heat are proportional to the density. As a consequence the cloud remains 
isothermal to a fairly good accuracy. 

Fig. [7] shows the change in the entropy per particle during the evolution of the system. 
The data points show the result of a simulation using an ideal gas equation of state with 
Eq/Ef = 1 and a shear viscosity fj = a n n with a n = 0.1. For the ideal gas equation of state 
we can compute the entropy using the Sackur- Tetrode formula. The dashed line shows the 
entropy absorbed by the the heat sink for a calculations with no reheating of the gas. In 
this case the produced entropy scales asymptotically as js] 



The full simulation tracks the results without reheating very well for (u>±t) < 3. At later 
times the full simulation produces less entropy then the model without reheating. However, 
we still find that the total entropy continues to grow as t — > oo. Numerically, we find that 
the asymptotic behavior is well described by (AS)/N ~ (w^t) 1 / 6 . 




(31) 
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FIG. 9: This figure shows the ratio of the transverse scale factor in dissipative hydrodynamics 
over the result in ideal hydrodynamics. The data points are from a simulation with r] ~ n 2 and 
(a n ) = 0.1. The lines show the same calculations as in Fig. |U 

C. Dependence on the equation of state and the functional form of the shear 
viscosity 

There are several arguments that indicate that the evolution of the system is not very 
sensitive to the equation of state P(n, T) as long as the universal relation P = |£ is satisfied. 
In Sect. HV] we showed that the exact solution of Euler's equation is independent of the 
equation of state. We also showed that the equations of dissipative hydrodynamics given in 
equ. 0231211) are independent of the equation of state as long as the velocity field remains 
exactly linear. In this section we will study the dependence on the equation of state using 
numerical simulations of the complete hydrodynamic equations. We will compare the results 
obtained using the ideal gas equation of state P = nT and the equation of state described 
in Appendix [A] We consider a temperature T = 0.25Tp, close to the superfluid phase 
transition, where the deviation of the experimental equation of state from the ideal gas 
equation is largest. For the ideal gas equation of state we have (E /Ep) = 3(T/Tp) = 
0.6. We choose a n = 0.06 so that f3 = 0.066. For the experimental equation of state we 
find Eq/Ef = 0.785 and we set d n = 0.0785 to keep [3 fixed. We find that the effect of 
the equation of state on the change in t cr is smaller than the accuracy of our calculation, 
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FIG. 10: This figure shows the evolution of the temperature profile in viscous hydrodynamics. 
The data points show the temperature T{x, 0,0, t) determined in a numerical simulation with 
r] ~ n 2 and (a) = 0.1 at several different times i = 0, 1.65, 2.66. The lines are the result in ideal 
hydrodynamics . 



We have also studied the dependence of dissipative effects on the functional form of the 
shear viscosity. The approximate solutions discussed in Sections HV] and [V] suggest that 
dissipative effects depend only on the trap average (a n ), see equ. ( 120]) . In the following we 
will test this idea by comparing calculations with a n ~ const, corresponding to rj ~ n, and 
&n ~ n/(mT) 3 / 2 , which implies t] ~ n 2 /(mT) 3 ' 2 . We write r] = ^n 2 / (mT) 3 / 2 and fix i] 2 
from (a n ). For a Gaussian profile 



Fig.[9]shows the evolution of the transverse radius for (a n ) = 0.1. The lines are the same as in 
Fig.Hl We observe that the calculations with rj ~ n and rj ~ n 2 are very similar for (u±t) < 3. 
At later times the non-linear dependence of rj on n leads to some extra acceleration. Fig. [TUI 
shows the corresponding temperature profiles. We observe that for rj ~ n 2 reheating takes 
place predominantly near the center of the cloud. This leads to an outward temperature 
gradient which is the source of the additional acceleration. Quantitatively, the difference 



(5t cr (P ld )-5t cr (P e *))/t cr < 10 



3 




(32) 
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between the r\ ~ n and rj ~ n 2 is about 25%, 

St cr (rj ~ n) - St cr (rj ~ n 2 ) 
<5t cr (r? ~ n) 



0.264 



(33) 



D. Rotating solutions 

In this section we study the time evolution of a rotating cloud. We take the initial velocity 
profile to be of the form v = aV(xz) with a(0) = u rot = 0Auj z [20J]. We have checked that 
the results are unaffected by taking the initial profile to be of the form v = Qy x x. The 
reason is that for a strongly deformed cloud the initial momentum density pv is essentially 
the same for irrotational or rigid initial conditions. 

We determine the angle of the major axis and the aspect ratio of the cloud. The angle is 
related to the Gaussian radii by 

2(xz) 



tan(20) 



{z 2 ) - (x 2 ) 



(34) 



and the aspect ratio is given by 



.4 



R 



(x 2 ) + (z 2 )+ ((z 2 ) - (x 2 )f + A{xz) 



1/2 >, 1/2 



1/2 



(35) 



(x 2 ) + (z 2 ) - [((z 2 ) - (x 2 )Y + A(xz) 2 

Our results are shown in Fig. [TTJ The solid line is the result in ideal hydrodynamics. A 
good approximation to the evolution of the angle in ideal hydrodynamics is 

a\ 2 b 2 b 2 



where b±,b z are the scale parameters for the pure expansion (without rotation) and 



(36) 



a(t) 



2u) ro tt 

A 2 



uj, < 1 , 



(37) 



•5»s£ wjl*>1, 

with 7 = 2/3. This result shows that the angle goes through 45 degrees at the same time 
at which the expanding system reaches an aspect ratio of 1. 

Fig. [TT1 shows the time evolution of the angle and the aspect ratio. The solid line shows the 
result in ideal hydrodynamics, the dashed line shows the result in dissipative hydrodynamics 
neglecting reheating, and the data points are from a numerical simulation with fj = d n n and 
a n = 0.1. We observe that the effect of reheating in rotating clouds is similar to the effect in 
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FIG. 11: This figure shows the time evolution of the angle (top panel) and the aspect ratio (bottom 
panel) of a rotating cloud as a function of time. The initial density profile is Gaussian, and the 
initial flow profile is an irrotational flow with co ro t = 0.4w 2 . The thin line is the result in ideal 
hydrodynamics, the dashed line is the result in dissipative hydrodynamics with a heat sink, and 
the thick line is result of a numerical calculation with a = 0.1. 

non-rotating systems. Reheating accelerates the system and reduces dissipative corrections. 
This effect can be quantified in terms of the time £45° at which the angle of the major 
axis passes through 45° (angular momentum conservation combined with the approximately 
irrotational nature of the flow implies that the aspect ratio never reaches the value 1). We 
find that, within the accuracy of our calculation, the dissipative correction to £45° is equal 
to the dissipative correction to the crossing time (see Sect. IVI Bj) . 5t^o = 5t cr . This implies, 



21 



in particular, that earlier estimates of the shear viscosity based on calculations that do not 
take into account reheating have to be a corrected by a factor ~ 2 Jg[ [22 ] . 



VII. CONCLUSIONS AND OUTLOOK 

In this work we studied the expansion dynamics of a dilute Fermi gas at unitarity in 
the framework of dissipative hydrodynamics. Our main goal was to study whether one 
can extract the shear viscosity from expanding systems. This is not immediately obvious, 
because in an expanding system all internal energy is eventually converted into kinetic 
energy, irrespective of whether there is dissipation or not. 

We find that shear viscosity does lead to characteristic effects in the expansion dynamics. 
Shear viscosity causes a characteristic curvature in the time evolution of the aspect ratio 
An(t) of the cloud. In ideal hydrodynamics internal energy is converted to kinetic energy 
very quickly, over a time period (u±t) < 3. After this time Au(t) is essentially linear. In dis- 
sipative hydrodynamics energy stored in the transverse motion is converted into longitudinal 
kinetic energy, and the longitudinal expansion takes place on a much longer time scale. As a 
result An(t) exhibits a characteristic curvature at times as large as (u±t) ~ A -1 ~ 25. This 
effect was recently observed by Cao et al. [23], which shows that dissipative hydrodynamics 
is indeed valid at (00 ±t) ~ 25. This is a remarkable discovery, because during the evolution 
the density drops by a factor A -2 ~ 10 3 . 

We also find that a quantitative description of the dependence of Ar(£) and other ob- 
servables on the shear viscosity has to include reheating. For a cloud with an aspect ratio 
of 25 the extracted shear viscosity is about a factor 2 too small if reheating is neglectec . 



This affects the estimates presented in 



22] but not the recent work of Cao et al. 23]. 



Reheating also does not affect estimates of the shear viscosity based on the damping of 



collective modes 
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19]. 



We showed that a determination of the shear viscosity does not require an accurate 
knowledge of the equation of state P(n,T). The only important aspect of the equation of 
state is the universal relation P = |£ . We also studied the dependence of viscous effects on 
the functional form of rj(n,T). We find that to first approximation the expansion dynamics 
constrains the cloud average of the shear viscosity. In this approximation the universal 
function a n (mT/n 2 ^ 3 ) can be determined by extracting (a n ) as a function of T/Tp from 
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data, and then inverting equ. (J21]). This only requires knowledge of the initial density 
profile. The result can be used as input for a more accurate determination based on full 
hydrodynamics. 

There are several issues that remain to be studied. The most important problem has 
to do with the breakdown of hydrodynamics in the dilute corona of the cloud. In the low 
density, high temperature limit the shear viscosity can be reliably computed. The result 
shows that the shear viscosity is independent of the density, 77 ~ (mT) 3 ' 2 . This implies that 
the total amount of heat dissipated by the dilute tail of the density distribution is infinite. 
We have previously argued that this problem can be resolved by taking into account the fact 
that the dissipative contribution to the stress tensor relaxes to the Navier-Stokes form on a 

n n 

time scale which is proportional to the density of the system [8|, |29[ . In kinetic theory we 
expect that TRdt(SHij) = (rjaij — SHij) where the relaxation time is given by tr = r)/(nT). 
This implies that in the dense regime the shear viscosity relaxes to its equilibrium value on 
a time scale that is fast compared to the time scale of the hydrodynamic expansion, but in 
the dilute regime dissipation is governed by an effective viscosity which is proportional to 
the density. 

This idea can be implemented by using an effective (a n ) in solving the equations of 
dissipative hydrodynamics js), l^j]. It is clearly preferable, however, to include the effects 
of finite relaxation time by including higher derivative terms in the equations of dissipative 
fluid dynamics (this is known as 2nd order, or Burnett, hydrodynamics), or by coupling the 
hydrodynamic description to kinetic theory. 
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Appendix A: Equation of state 



In this appendix we describe a parameterization of the e quat ion of state in the normal 



phase. The equation o 



Carlo simulations 
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of Nascimbene et al. 



321 ] , using quantum Monte 



: state has been studied experimentally 130 

nn " - 



35], and many-body theory [36|, |37|. Here we follow the recent work 



3l| and write 



(Al) 
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FIG. 12: Equation of state of the unitary Fermi gas in the normal phase. In this figure we show 
the function h((), where £ is the fugacity. The data points are from Nascimbene et al., the thick 
solid line shows the parameterization discussed in the text, the dashed line is the non- interacting 
gas result h = 1, and the dotted line shows the second order Virial expansion. 

where Pi (/i, T) is the ideal gas equation of state of a single species non-relativistic Fermi gas 

P l {^T) = -T\f B Li, /2 {-C l ), (A2) 

and XdB = [(27r)/(mT)] 1 / 2 is the de Broglie wave length. Here, Li a (x) is the Polylogarithm 
function, and ( = exp(— fi/T) is the fugacity. We parameterize h(() as 

ho e + cic + c 2 



C 2 + C 3 C + C 4 



(A3) 



and determine the parameters q from a fit to the data of Nascimbene et al. [3j|. This 
parameterization is motivated by the fact that the data for ( > 1 is very well described 
by the Virial expansion h(()/2 = 1 + b 2 /( + h/( 2 + 0(l/( 3 ). At unitarity b 2 = l/y/2 and 
6 3 = | — 0.355 = 0.23. The value of h(() at zero fugacity is related to the Bertsch parameter 
£ = ^/E F . Using £ ~ 0.4 we have h(0)/2 = £~ 3 / 2 ~ 3.8. A fit for fugacities in the range 
C G [0.03, 5] gives 

ci = 1.3543, c 2 = -0.0174, c 3 = 0.5724, c 4 = -0.0084. (A4) 

We compare the data to our fit and the Virial expansion in Fig. [12j From the pressure we 
can determine other thermodynamic quantities. The density and entropy density are given 
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FIG. 13: Density of a trapped Fermi gas in the local density approximation at a temperature 
T/Tp = 0.25, slightly above the critical temperature. The thick line is the result based on the 
equation of state from Nascimbene et al., the dashed line is the result for a free gas, and the dotted 
line shows the result using the equation of state in the high temperature limit. 

by 

n(jM, T) = X d lg(C) , s(jm, T) = \-*k({) , (A5) 

with 

g(0 = -Lh /2 {-C l )h{C) + CLkfzi-C^h'iC) , (A6) 
HO = ~ (log(0^s/a(-C" 1 ) + ^B/aC-C" 1 )) HO 

+ \og(C)Lz 5/2 (-C 1 )h'(0- (A7) 

In a trapped system we use the local density approximation n(x) = n(/i(x),T) with fj,(x) = 
fi — V(x) where V(x) is the trapping potential. This determines the density profile if the 
temperature and the chemical potential (or the fugacity) at the center of the trap are given. 
In practice we usually specify the temperature and the total number of particles. The 
particle number defines a temperature scale Tp = (3iV) 1//3 a), where u = {uo x ujyUj z ) 1 ^ is the 
geometric mean of the trap frequencies. Given T/Tp the fugacity Co at the center of the 
trap is determined by the condition 

P^(^)7^(<o«p(t))-i- (as) 

This equation has to be solved numerically. In the high temperature limit Co = 6(T/Tp) 3 . 
In Fig. 1151 we show the density profile at T/T F = 0.25. We show the exact density, the 
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density of a free gas, and the high temperature (Gaussian) approximation. The effects of 
quantum degeneracy decrease the central density, whereas interactions increase the density. 
The two effects partially cancel and the exact density is about 50% larger than the Gaussian 
approximation. 

Once the initial density and pressure have been determined the equations of fluid dy- 
namics fix the evolution of P and n. The equation of state is needed in order to compute 



other thermodynamic quantities like the temperature and the chemical potential 38[. The 
fugacity can be computed from 

2 f^\ V2 _ 2 /(0 3/2 = F(n (m 

(2tt) 3 / 2 U 5 /V K.f(C)] 5/2 

where /(C) = -Li h/2 {-C,~ l )h{C). E Q U - O implies that 



C 



I 2 m 3/2p3/2 \ 

1^)3/2 n 5/2 ) ■ ( A1Q ) 



In general, F 1 (y) has to computed numerically. In the high temperature limit F 1 (y) ~ y. 
Once the fugacity is known the temperature can be computed from 

Cf'(C) P 

T = --J77V-- ( A11 ) 
/(C) n 

In the high temperature limit /(C) — 2/C which implies C/'(C)//(C) — — 1 an d T = P/n. 
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